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ABSTRACT 

We study the mass, velocity dispersion, and anisotropy profiles of ACDM halos using 
a suite of N-body simulations of unprecedented numerical resolution. The Aquarius 
Project follows the formation of 6 different galaxy-sized halos simulated several times 
at varying numerical resolution, allowing numerical convergence to be assessed directly. 
The highest resolution simulation represents a single dark matter halo using 4.4 billion 
particles, of which 1.1 billion end up within the virial radius. Our analysis confirms a 
number of results claimed by earlier work, and clarifies a few issues where conflicting 
claims may be found in the recent literature. The mass profile of ACDM halos deviates 
slightly but systematically from the form proposed by Navarro, Frenk & White. The 
spherically-averaged density profile becomes progressively shallower inwards and, at 
the innermost resolved radius, the logarithmic slope is 7 = — dlnp/dlnr £ 1. Asymp- 
totic inner slopes as steep as the recently claimed p oc r~ 12 are clearly ruled out. 
The radial dependence of 7 is well approximated by a power-law, 7 oc r a (the Einasto 
profile). The shape parameter, a, varies slightly but significantly from halo to halo, 
implying that the mass profiles of ACDM halos are not strictly universal: different 
halos cannot, in general, be rescaled to look identical. Departures from similarity are 
also seen in velocity dispersion profiles and correlate with those in density profiles 
so as to preserve a power-law form for the spherically averaged pseudo-phase-space 
density, p/a 3 oc r -1 - 8 ' 5 . The index here is identical to that of Bertschinger's sim- 
ilarity solution for self-similar infall onto a point mass from an otherwise uniform 
Einstein-de Sitter Universe. The origin of this striking behaviour is unclear, but its 
robustness suggests that it reflects a fundamental structural property of ACDM halos. 
Our conclusions are reliable down to radii below 0.4% of the virial radius, providing 
well-defined predictions for halo structure when baryonic effects are neglected, and 
thus an instructive theoretical template against which the modifications induced by 
the baryonic components of real galaxies can be judged. 
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1 INTRODUCTION 

A couple of decades of steady progress in the simulation of 
non-linear structures in a cold dark matter (CDM) domi- 
nated universe have resulted in significant advances in our 
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understanding of the clustering of dark matter on the scale 
of galactic halos. There is now widespread consensus that 
the hierarchical assembly of CDM halos yields: (1) mass pro- 
files that are approximately "universal" (i.e., independent of 
mass and cosmological parameters aside from simple physi- 
cal scalings (Navarro et al., 1996, 1997, hereafter NFW), (2) 
strongly triaxial shapes, with a slight preference for nearly 
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prolate systems (e.g., Frenk et al., 1988; Jing & Suto, 2002; 
Allgood et al., 2006; Hayashi et al., 2007), (3) abundant, 
but non-dominant, substructure (Klypin et al., 1999; Moore 
et al., 1999a; Ghigna et al., 2000; Gao et al., 2004), and 
(4) "cuspy" inner mass profiles, where the central density 
increases systematically as the numerical resolution of the 
calculation is improved (see, e.g., NFW, Moore et al., 1999b; 
Fukushige & Makino, 2001; Navarro et al., 2004; Diemand 
et al., 2005). 

Despite this consensus, there are a number of issues 
where conflicting claims may be found in the recent liter- 
ature, hindering the design and interpretation of observa- 
tional tests aimed at validating or ruling out various aspects 
of the CDM theory on these scales. One contentious issue 
concerns the statistics, spatial distribution, and structure of 
substructure, and their consequences for the discovery and 
interpretation of possible signals of dark matter annihilation 
in the gamma-ray sky (Stoehr et al., 2003; Diemand et al., 
2007; Kuhlen et al., 2008; Springel et al., 2008a, and refer- 
ences therein). The controversy extends to the structure of 
the inner cusps both of the main halo and of substructure 
halos, where some recent work has claimed a well-defined 
central slope of p oc r -1 ' 2 (Diemand et al., 2004, 2005, 2008) 
whereas others have argued that no compelling evidence for 
such power-law behaviour is apparent (Navarro et al., 2004; 
Graham et al., 2006). 

Considerable debate also surrounds whether the struc- 
ture of CDM halos is truly "universal". This is indeed the 
case if halos have mass profiles that are well-described by 
two-parameter formulae, such as the NFW profile or some 
of its modifications; see, for example, Moore et al. (1999b, 
hereafter M99). These profiles have two scaling parameters 
(mass and size) but fixed shape, so that two different halos 
can, in principle, be rescaled to be indistinguishable from 
each other. 

On the other hand, recent work suggests that at least 
three parameters may be needed to describe halo mass pro- 
files accurately. An example is the Einasto formula (Einasto, 
1965), shown by Navarro et al. (2004) to improve signif- 
icantly the accuracy of the fits to the inner density pro- 
files of simulated halos. It is unclear from that work, how- 
ever, whether the improvement is due to the fact that the 
Einasto formula has a different asymptotic inner behaviour 
than NFW or to the extra shape parameter it introduces. 
Merritt et al. (2005, 2006) explored this further and argued 
that the third parameter is indeed needed to account faith- 
fully for the curvature in the shape of the density profile. 
Merritt et al.'s conclusions have received support from the 
work of Gao et al. (2008) and Hayashi & White (2008), who 
have stacked density profiles of many halos of similar mass to 
show that mean profile shape, and, in particular, the Einasto 
shape parameter a (see eq. 4 below), depends systematically 
on halo mass. This implies that the mass profile of ACDM 
halos is not strictly universal; no simple scaling of the av- 
erage profile of cluster halos will provide an accurate fit to 
the average profile of galaxy halos. 

Many of these controversies and uncertainties may be 
traced to the fact that earlier work has lacked the numerical 
resolution and the representative halo sample needed to set- 
tle the debate. For example, the dark matter annihilation 
flux observable from Earth depends crucially on resolving 
not only substructures but also the nested "substructure 



within substructure" expected from the hierarchical assem- 
bly of CDM halos. Only the most recent simulations have 
been able to begin addressing this issue (see, e.g., Diemand 
et al., 2008; Springel et al., 2008a,b). 

A similar comment applies to the structure of the inner 
cusp, where pinning down the asymptotic inner behaviour of 
the dark matter density profile depends crucially on under- 
standing the limitations introduced by, for example, finite 
particle number, gravitational softening, and time-stepping 
technique. 

We have shown in earlier work (Power et al., 2003, here- 
after P03) that, when suitable choices of the numerical pa- 
rameters are made, the main factor determining the inner- 
most radius where the mass profile may be measured reli- 
ably is the total number of particles used in the simulation. 
Empirically, the boundary of the region where numerical 
convergence is achieved roughly corresponds to the radius 
where the two-body relaxation time, irciax, exceeds the age 
of the Universe. Since i rc iax scales roughly like the enclosed 
number of particles times the local orbital timescale, and 
the latter drops sharply toward the centre, extending the 
resolved region inwards even modestly requires a dramatic 
increase in the total number of particles. 

These difficulties, coupled to the significant halo-to-halo 
scatter already seen in early work, imply that substantive 
progress on these issues requires a concerted numerical ef- 
fort where several different halos are simulated with varying 
numerical resolution, so that cosmic variance and numerical 
convergence may be assessed directly. 

These are the aims of The Aquarius Project, a recently 
completed suite of numerical simulations of the formation of 
galaxy-sized halos in the ACDM cosmogony. The series in- 
cludes re-simulations of six different ~ 10 12 Mq halos where 
the number of particles is systematically varied. In one case, 
the same halo is simulated 5 times, increasing gradually the 
number of particles in the halo from about one million to 
~ 1.1 billion within the virial radius. The highest resolution 
simulations of the other 5 halos have roughly 100-200 million 
particles each within the virialized region. 

The simulation series has been presented recently by 
Springel et al. (2008a, b), where the interested reader may 
find relevant details. Our first paper (Springel et al., 2008a) 
deals with predictions of the annihilation signal whereas the 
second (Springel et al., 2008b) addresses the statistics, spa- 
tial distribution, and structure of dark matter substructures. 
Here we deal with the structure of the main halo, with spe- 
cial emphasis on the structure of the inner cusp. The plan 
of the paper is as follows. Sec. 2 summarizes briefly the nu- 
merical parameters of our simulations; Sec. 3 and 4 present 
our main results. We conclude with a brief discussion and 
summary in Sec. 5. 



2 THE NUMERICAL SIMULATIONS 

We present here for completeness a brief summary of the 
numerical simulations, and refer the reader to Springel et 
al. (2008a,b) for further details. 
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Figure 1. Spherically-averaged density (left) and circular velocity (right) profiles for the Aq-A halo simulation series. Different colours 
correspond to different resolution runs, as labeled in the figure. The density profile is multiplied by r 2 in order to emphasize small 
deviations. The bumps in the outer regions may be traced to the presence of substructure and unrelaxed tidal debris. Profiles are shown 
from ~ 3r2oo down to the "convergence radius", r'onvi corresponding to the radius where the relaxation time, t ro i ax , is of the order of 



(7) - 

the age of the Universe. The thick portion of each profile indicates the region r > r^ a ' nv where t re iax is more than 7 times the age of the 

(7) 

universe and where stricter convergence is achieved. Outside r^ onv circular velocity estimates converge to better than 2.5% (see Fig. 2). 
The dot-dashed line shows an Einasto profile with a = 0.17 matched at (r—2,p— 2), the peak in the r 2 p profile. This provides an excellent 
fit to the structure of the inner regions of the halo, as shown by the residuals plotted in the bottom panels. Arrows indicate the softening 
length h s of each simulation. 



2.1 The Cosmological Parameters 

All our simulations assume a ACDM cosmogony with the 
following parameters: fl m = 0.25, £Ia = 0.75, as = 0.9, 
n s — 1, and Hubble constant Ho = 100 h kms -1 Mpc -1 = 
73 kms^ 1 Mpc -1 . These cosmological parameters are the 
same adopted in previous numerical work by our group, 
such as the Millennium Simulation of Springel et al. (2005) , 
and are consistent, within their uncertainties, with con- 
straints derived from the WMAP 1- and 5-year data analy- 
ses (Spergel et al., 2003; Komatsu et al., 2008) and with the 
recent cluster abundance analysis of Henry et al. (2008) . 

2.2 The Code 

The simulations were carried out with a new version of 
the GADGET (Springel et al., 2001; Springel, 2005) parallel 
cosmological code. This version, which we call GADGET-3, 
has been especially developed for this project, and imple- 
ments a novel domain decomposition technique in order to 
achieve unprecedented dynamic range in massively-parallel 
computer systems without sacrificing load balancing or nu- 
merical accuracy. Time stepping is carried out with a kick- 
drift-kick leap-frog integrator where the timesteps are based 
on the local gravitational acceleration, together with a con- 
servatively chosen maximum timestep for all particles. 

Pairwise particle interactions are softened with a spline 
of scalelength h s , so that they are strictly Newtonian for par- 



ticles separated by more than h s . The resulting softening is 
roughly equivalent to a traditional Plummer-softening with 
scalelength ec ~ h s /2.8. The gravitational softening length 
is kept fixed in comoving coordinates throughout the evo- 
lution of all our halos. The dynamics is then governed by a 
Hamiltonian and the phase-space density of the discretized 
particle system should be strictly conserved as a function of 
time (Springel, 2005). 



2.3 Halo Selection 

All halos in the Aquarius suite were identified for resimu- 
lation in a 900 3 -particle parent simulation of a 100 /i -1 Mpc 
box. The identification technique selects all ~ 10 12 M© halos 
in the box and chooses, at random, a few of them that sat- 
isfy a mild isolation criterion (no neighbour exceeding half 
its mass within 1ft -1 Mpc). This criterion is only imposed 
in order to remove halos in the vicinity of massive groups 
and clusters, which may have evolved differently from the 
average. 

Each halo is then resimulated at various resolutions, 
making sure that each resimulation shares the same power 
spectrum and phases at all resolved spatial frequencies. Ini- 
tial displacements are imprinted using the Zeldovich ap- 
proximation and a 'glass-like' uniform particle load (White, 
1996). The 100ft _1 Mpc simulation box is divided into 
a "high-resolution" region, which corresponds to the La- 
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Figure 2. Top panel: Fractional deviations in the circular velocity 
profile of the Aq-A convergence series versus the (enclosed) relax- 
ation time, tj-eiaxi expressed in units of the circular orbit period at 
the virial radius, t irc(f20o)' Deviations are measured relative to 
the highest resolution halo, Aq-A-1. Note that departures from 
convergence for all simulations are similar when expressed this 
way, indicating that t re iax is the main parameter determining 
convergence. Solid circles mark the location of the convergence 
criterion proposed by P03. Note that V c estimates converge there 
to about 10%. A stricter convergence criterion, e.g., 2.5% conver- 
gence in V c , is achieved at larger radii, where t re iax ~ 7t c i rc (r2oo) 
(right vertical line). Bottom panel: Relaxation time versus radius 
for all five Aq-A simulations. Arrows indicate h s = 2.8 eg, the 
lengthscale where pairwise interactions become Newtonian. 



grangian region surrounding the target halo, and a low- 
resolution region (the rest of the box), which is represented 
with a smaller number of particles with mass increasing with 
distance to the target halo. We have carefully designed the 
geometry of the high-resolution region in order to avoid con- 
tamination of the halo by massive low-resolution particles. 
Typically, about 30% of particles in the high-resolution re- 
gion end up in the virialized region of the final halo, and no 
higher mass particles end up within the virial radius of the 
final halo. 

Table 1 lists some basic information about each simu- 
lation. This includes a symbolic simulation name, the parti- 
cle mass in the high-resolution region, m p , the gravitational 
softening length, ec the virial radius*, r2oo, as well as the 
total mass, M200 and the total number of particles, N200, 
enclosed within r2oo- Other structural parameters of inter- 
est include the location of the peak in the circular velocity 
profile, specified by r max and Vm ax , as well as that of the ve- 
locity dispersion profile (<r max and r(er max )). <7host indicates 
the ID rms velocity of the main halo within r2oo (excluding 
substructures). 

Table 1 lists only information on the halos used in this 
paper. A more complete list of numerical parameters may 
be found in Springel et al. (2008b). One of our halos, la- 
beled Aq-A, has been resimulated 5 times, spanning a fac- 
tor of ~ 2000 in particle mass. Our naming convention uses 
the tags "Aq-A" through "Aq-F" to refer to each of the six 
Aquarius halos. An additional suffix "1" to "5" denotes the 
resolution level. "Aq-A-1" is our highest resolution calcula- 
tion: it follows the surroundings of Aq-A with ~ 4.4 billion 
particles, ~ 1.1 billion of which end up within T2oo- We have 
level-2 simulations of all 6 halos, corresponding to between 
100 and 200 million particles per halo (within r2oo)- The 
softening parameters of each simulation adopt the "optimal" 
softening recommendation of P03, which aims to balance the 
number of timesteps required for accurate integration whilst 
minimizing the loss of spatial resolution. 



2.4 Radial Profiles 

Our analysis uses spherically-averaged profiles of the ba- 
sic dynamical properties describing the structure of ACDM 
halos: the density, circular velocity, velocity dispersion, 
and anisotropy profiles. Typically, these are computed in 
50 spherical shells equally spaced in log 10 r (where r is 
the distance to the halo center), and spanning the range 

1.5 x 10~ 4 < r/r2oo < 3. (When different choices for ei- 
ther the number of bins or the radial range are made, this 
is stated explicitly in the analysis below.) These concentric 



* We define the virial mass of a halo, M200, as that contained 
within a sphere of mean density 200 X p cr it. The virial mass de- 
fines implicitly the virial radius, r200i and virial velocity, V200 = 
(GAf20o/^20o) 1//2 , of a halo, respectively. We note that other defi- 
nitions of "virial radius" have been used in the literature; the most 
popular of the alternatives adopts a density contrast (relative to 
critical) of A rj 178 H^ 45 ~ 100 (for our adopted cosmological 
parameters, see Eke et al. (1996)). We shall refer to these alter- 
native choices, where appropriate, with a subscript indicating the 
value of A; i.e., r-50 would be the virial radius obtained assuming 
A = 50, and so an enclosed density 200 times the mean cosmic 
value. 
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Figure 3. Left: Spherically-averaged density profiles of all level-2 Aquarius halos. Density estimates have been multiplied by r 2 in order 
to emphasize details in the comparison. Radii have been scaled to r_2, the radius where the logarithmic slope has the "isothermal" value, 
—2. Thick lines show the profiles from r'^, outward; thin lines extend inward to rconv For comparison, we also show the NFW and M99 
profiles, which are fixed in these scaled units. This scaling makes clear that the inner profiles curve inward more gradually than NFW, 
and are substantially shallower than predicted by M99. The bottom panels show residuals from the best fits (i.e., with the radial scaling 
free) to the profiles using various fitting formulae (Sec. 3.2). Note that the Einasto formula fits all profiles well, especially in the inner 
regions. The shape parameter, a, varies significantly from halo to halo, indicating that the profiles are not strictly self-similar: no simple 
physical rescaling can match one halo onto another. The NFW formula is also able to reproduce the inner profiles quite well, although 
the slight mismatch in profile shapes leads to deviations that increase inward and are maximal at the innermost resolved point. The 
steeply-cusped Moore profile gives the poorest fits. Right: Same as left, but for the circular velocity profiles, scaled to match the peak of 
each profile. This cumulative measure removes the bumps and wiggles induced by substructures and confirms the lack of self-similarity 
apparent in the left panel. 



shells are centered at the location of the particle identified by 
the SUBFIND algorithm (Springel et al., 2001) as having the 
minimum gravitational potential. Extensive tests show that 
this procedure identifies the region where the local density of 
the main subsystem of each halo peaks, and is coincident in 
most cases (except perhaps major ongoing mergers between 
comparable-mass halos) with the results of other methods, 
such as the "shrinking sphere" method discussed by P03. 



The mass density in each radial bin is estimated as the 
dark mass in the bin divided by its volume, and assigned to 
a radius corresponding to the bin center. Circular velocities 
are computed by adding up the mass of each bin plus all 
interior ones, and assigned to the radius corresponding to the 
outer edge of the bin. The construction of velocity dispersion 
and anisotropy profiles is described in detail in Sec. 4.1. 
When differentiation is necessary, such as when computing 
the logarithmic slopes shown in Figs. 5 and 6, we use a simple 
3-point Lagrangian interpolation to perform the numerical 
differentiation (as implemented by the DERIV subroutine of 
the IDL software package). 



3 MASS PROFILES 

3.1 Numerical Convergence 

We begin our study of the mass profile by using our series 
of re-simulations of the Aq-A halo in order to assess the ra- 
dial range where numerical convergence is achieved. Figure 1 
shows the mass profile of the five Aq-A resimulations; the left 
panels show the spherically-averaged density profile (multi- 
plied by r 2 in order to emphasize small departures) ; the right 
panels the corresponding circular velocity profile. Lines of 
different colours correspond to different resimulations, as la- 
beled. Arrows indicate h s = 2.8 ec, the lengthscale where 
softened pairwise interactions become fully Newtonian. 

This figure demonstrates the striking numerical conver- 
gence achieved in our re- simulations. Outside some char- 
acteristic radius (which we discuss below), all the profiles 
are essentially indistinguishable from each other, even down 
to details such as "bumps" in the outer regions caused by 
the presence of substructure. As discussed by Springel et al. 
(2008b), this reflects the high quality of the numerical in- 
tegration of GADGET-3 and the careful approach we have 
taken to building our initial conditions; indeed, the Aq-A 
resimulations not only reproduce faithfully the properties 
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of the main halo, but even the mass, location and internal 
structure of most major substructures. 

Inevitably, near the centre the mass profiles diverge as a 
consequence of numerical limitations. Each profile is plotted 
down to the "convergence radius" proposed by P03. These 
authors demonstrate that deviations from convergence de- 
pend (for appropriate choices of other numerical parame- 
ters) solely on the number of particles, and scale roughly 
with the collisional "relaxation" time, i rc iax- Expressed in 
units of the circular orbit timescale at r2oo (which is of the 
order of the age of the Universe), k = Wax Acirc^oo), the 
relaxation time may be written as: 

<r) _ N r/V c _ y/200 ^ N(r) _ f g(r) V 1/2 ^ (1) 



8 In N r2oo /V200 



8 lniV(r) 



where N = N(r) is the enclosed number of particles and 
p(r) is the mean enclosed density within r. 

According to P03, deviations of roughly 10% are ex- 
pected in the V c profile where and they adopted this 
condition to define the convergence radius, r conv . Stricter 
convergence demands larger values of k, and we shall use 
a superscript on r conv to denote the value of k adopted for 
its definition. For instance, rconv 

K = 1. 



»conv corresponds to 



Profiles in Fig. 1 are thus plotted in the range 



3r2oo]- As shown in the bottom right panel, this inner radius 
indeed corresponds to the point where systematic deviations 
in Vc(r) reach ~ 10%. It is also clear from this figure that 
convergence in the local density profile is always much easier 
to achieve, so concentrating our analysis on the enclosed 
mass profile, or on the circular velocity, is a conservative 
approach. 

Although each halo converges over a different radial 
range, the departures from convergence are all similar when 
expressed in terms of n. This is shown in the top panel in 
Fig. 2, where differences in V c from our highest-resolution 
halo, Aq-A-1, are shown as a function of k, for the other Aq- 
A resimulations. Deviations of ~ 10% are typical at k = 1; 
convergence to better than ~ 2.5%, on the other hand, re- 
quires k ~ 7 (indicated by the right dashed vertical line). 

We may use these results to estimate convergence radii 
for our highest resolution run, Aq-A-1: its V c profile con- 
verges to better than 10% for radii r > r conv = 112 /i -1 
pc; 2.5% convergence or better is expected for r > rSiv = 
253 pc (see bottom panel of Fig. 2). Convergence radii 
for various values of n are listed in Table 2 for each simulated 
halo. 



3.2 Fitting formulae 

The fitting formulae we have used to describe the mass pro- 
file of our simulated halos are the following: (i) The NFW 
profile, given by 

Ps 



p(r 



(2) 



(r/r s )(l+r/r s ) 2 ' 
(ii) the modification to the NFW profile proposed by M99 
pM 



P{r) {r/r M y- s [l + {r/r M Y' 5 ] 
and (iii) the Einasto profile, 

In(p(r)/p_ 2 ) = (-2/a)[(r/r_ a )° 



(3) 
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Figure 4. Minimum-Q values as a function of the Einasto param- 
eter a for best fits to all level-2 halo profiles in the radial range 
0.01 < r/r_2 < 5. Colors identify different halos, and line types 
the number of bins chosen for the profile. The minimum-Q val- 
ues obtained for NFW and M99 best fits are also shown, and arc 
plotted at arbitrary values of a for clarity. Note that Einasto fits 
are consistently better than NFW which are consistently better 
than M99, and that a significant improvement in Q is obtained 
when letting a vary in the Einasto formula. Q is approximately 
independent of the number of bins used in the profile, and is min- 
imized for different values of a for each individual halo. See text 
for further details. 



Because each of these formulae defines the character- 
istic parameters in a slightly different way, we choose to 
reparametrise them in terms of r_2 and p_2 = p(r_2), which 
identify the "peak" of the r 2 p profile shown in the left panel 
of Fig. 1. This marks the radius where the logarithmic slope 
of the profile, 7(1*) = — d In p/dlnr, equals the isothermal 
value, 7 = 2. 

The characteristic radius, r_2, is a well-defined scale- 
length which is relatively easy to identify in each halo with- 
out resorting to any particular fitting formula. In practice, 
we determine r_2 by computing the logarithmic slope pro- 
file, y(r), and identifying where a low-order polynomial fit 
to it intersects the isothermal value. Each r 2 p profile is then 
visually inspected in order to ensure that r_2 corresponds 
to the main peak of the profile, and that it is not unduly 
influenced by secondary peaks that arise as the result of sub- 
structure. (See the left panel of Fig. 1.) Table 2 lists r_2 and 
p-2 for all our simulated halos. Note that for the NFW pro- 
file, r_2 = r s and p_2 = Ps/4, while for the Moore profile, 
p-2 = (4/3) p M and r_ 2 = 2~ 2/3 r M - 

We note that, unlike NFW or M99, when a is al- 
lowed to vary freely the Einasto profile is a 3-parameter 
fitting formula. This is not, of course, the only possible 
extension of NFW-like profiles which allows for a variable 
shape with the aid of an extra free parameter. For example, 
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Figure 5. Logarithmic slope of the density profile as a function 
of radius for our Aq-A convergence series. As in other plots, thick 
lines show results for r > r^ Q ' nv , thin lines extend the profiles down 
to the less strict convergence radius r^j, v . Comparison shows that 
excellent numerical convergence for the slope is achieved down 
to a radius intermediate between these two convergence radii. 
Applied to the highest-resolution Aq-A-1 simulation, this implies 
that the slope is shallower than the asymptotic value of the NFW 
profile (r _1 ) in the inner regions. We see no sign of convergence to 
an asymptotic inner power-law. Instead, the profiles get shallower 
toward the centre as predicted by the Einasto formula (a straight 
line in this plot). The "critical solution" of Taylor & Navarro 
(2001) (which has a r~ °- 75 asymptotic inner cusp) does better 
than NFW but not as well as Einasto in reproducing the inner 
profile of the halo. 



Merritt et al. (2006) compared N-body halos with the 3- 
parameter Einasto formula, as well as with the anisotropic 
model of Dehnen & McLaughlin (2005) and with the de- 
projected Sersic (1968) model of Prugniel & Simien (1997). 
Merritt et al. conclude that, overall, Einasto's formula per- 
forms best. Therefore, we adopt it here for the rest of our 
analysis, although we do not exclude the possibility that 
other 3-parameter formulae may perform at least as well as 
Einasto's. A full exploration of this issue is beyond the scope 
of this paper. 



3.3 Fitting procedure 

Best-fit parameters are found by minimizing the deviation 
between model and simulation across all bins in a specified 
radial range. In the case of the density profile, the best fit 
is found by minimizing the figure-of-merit function, Q 2 , de- 
fined by 

JV bins 

Q 2 = J7— V(mpi-lnpr de1 ) 2 - (5) 

JVbins * ' 

i=l 
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log r/r_ a 



Figure 6. As Fig. 5, but for all level-2 resolution Aquarius halos, 
after scaling radii to r_2- 



This function provides an intuitively simple measure 
of the level of disagreement between simulated and model 
profiles. It is dimensionless; it weights different radii loga- 
rithmically; and, for given radial range, Q 2 is approximately 
independent of the number of bins used in the profile. Thus, 
minimizing Q 2 yields for each halo well-defined estimates of 
a model's best-fit parameters. Note that when Q is small, 
it is just the rms fractional deviation of the data from the 
model. 

It is less clear how to define a goodness-of-fit mea- 
sure associated with Q 2 and, consequently, how to assign 
statistically-meaningful confidence intervals to the best-fit 
parameter values. This difficulty arises because, at the very 
high resolution of the simulations analyzed here, discreteness 
noise in the binned density estimates is negligible. The figure 
of merit of a fit therefore depends not only on how faithfully 
a model approximates a halo but also on the presence of indi- 
vidual halo features that no simple fitting formula can hope 
to reproduce. These distinct features are present on small 
scales (substructure) and large scales (such as streams, as- 
phericity, and other relics of each halo's specific assembly 
history; see, e.g., Vogelsberger et al., 2009). As a result, bin- 
to-bin residuals are distinctly non-Gaussian and highly cor- 
related, precluding the use of simple statistical tools such as 
the v 2 distribution in order to assess goodness of fit. 

Assessing the acceptability of various Q values would 
require the definition of a detailed statistical model in order 
to measure reliably the departures of individual halos from 
a smooth profile whose average shape (and scatter) could be 
obtained directly by averaging various numerical realizations 
of halos of the same mass. Unfortunately, such procedure is 
unlikely to be robust with only 6 halos in our sample. 

Therefore, we limit our analysis to comparing the 
minimum-Q values obtained with various formulae, and to 
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discussing how Q changes as the fitting parameters are var- 
ied. The actual value of Q is, after all, a reliable and objec- 
tive measure of the average per-bin deviation from a par- 
ticular model. As we discuss below, this is, in many cases, 
enough to prefer unequivocally one fitting formula over an- 
other and to make a compelling case for the need of an extra 
parameter in the fit. 



3.4 Einasto vs NFW vs M99 

The left panel of Fig. 3 compares the density profiles of 
all six level-2 Aquarius halos, after scaling radii to r_2 and 
densities to p_2- The right-hand panel shows the circular 
velocity profiles, scaled in an analogous manner to match 
the peak of the profile, identified by r max and Knax- In these 
scaled units, the fitting formulae introduced in Sec. 3.2 are 
curves of fixed shape and normalization, as shown by the 
thin solid, dashed, and dot-dashed curves in Fig. 3. (The 
Einasto curve adopts a — 0.159 in this figure.) 

Comparison with the simulations (thick curves) indi- 
cates that there is a clear mismatch between the shape of 
the halo profiles and those of the NFW and M99 fitting 
formulae. This is not just a result of enforcing the r_2-p-2 
scaling. We illustrate this by showing, in the two bottom 
panels of Fig. 3, residuals from best fits obtained by adjust- 
ing both fit parameters of the NFW and M99 profiles (r_2 
and p-2) in order to minimize Q 2 . (The radial range chosen 
for these fits is r™ < r < 0.5r2oo-) Note the "S" shape 
in the residuals, which are largest (and increasing) at the 
innermost radius of the profile. Because of the shape mis- 
match, extrapolating either the NFW or M99 fits further 
inwards, to regions less well resolved numerically, is almost 
guaranteed to incur substantial error. 

The large-scale radial trend of the residuals from the 
best Einasto fits (middle panels of Fig. 3) , on the other hand, 
is rather weak, suggesting that the shape of the simulated 
halo profiles are much better accommodated by this formula. 
This is not just a result of the extra shape parameter in 
the Einasto formula: even when keeping a fixed to a single 
value, residuals are smaller and have less radial structure 
than those from either NFW or M99. 

We show this in Fig. 4, where we plot the minimum- 
Q (Qmin) values of the best Einasto fits for all six level-2 
Aquarius halos, as a function of the shape parameter a. For 
given value of a the remaining two free parameters of the 
Einasto formula are allowed to vary in order to minimize 
Q 2 . Different line types correspond to different numbers of 
bins used to construct the profile (from 20 to 50), chosen to 
span in all cases the same radial range, 0.01 < r/r_2 < 5, 
a factor of 500 in radius. Minimum- Q values are computed 
using a similar procedure for the NFW and M99 formulae, 
and are shown, for each halo, with symbols of corresponding 
colour. 

In terms of Q m i n , Einasto fits are consistently superior 
to NFW or M99, whether or not the a parameter is adjusted 
freely. For example, for fixed a = 0.15, all Einasto best fits 
have minimum-Q values below ~ 0.03. For comparison, best 
NFW and M99 fits have an average (Q min ) ~ 0.06 and 0.095, 
respectively. These numbers correspond to A^ms = 20, but 
they are rather insensitive to AWs, as may be judged from 



the small difference between the various lines corresponding 
to each halo in Fig. 4. 

We emphasize that, although the improvement obtained 
with Einasto's formula is significant, NFW fits are still ex- 
cellent, with a typical rms deviation of just ~ 6% over a 
range of 500 in radius. The use of the NFW formula may 
thus be justified for applications where this level of accuracy 
is sufficient over this radial range. 

When a is adjusted as a free parameter, (Qmin) ~ 0.018 
for Einasto fits. Furthermore, there is, for each halo, a well 
defined value of a that yields an absolute minimum in Q. 
The Q-dependence on a about this minimum is roughly sym- 
metric and, as expected, nearly independent of the number 
of bins used in the profile. The minimum in Q is sharp; 
a shift of just 0.015 in a typically leads to an increase of 
~ 50% in Q around the minimum. Given that the value of a 
that minimizes Q varies from 0.130 for Aq-E-2 to 0.173 for 
Aq-B-2, we conclude that the improvement obtained when 
allowing a to vary is significant. We quote nominal "error 
bars" for a in Table 2 that bracket the interval where Q 
deviates by less than 50% from the absolute minimum in 
Fig. 4. 

3.5 Self-similarity? 

The need for a variable a discussed above illustrates one 
of our main findings: namely, that the mass profiles of our 
Aquarius halos are not strictly self similar. The shapes of the 
profiles are subtly but significantly different from each other, 
and no rescaling can match one exactly to another. Halo 
Aq-E-2 provides the most striking example, deviating from 
halo Aq-D-2, for example, by almost a factor of 2 in density 
at ~ 0.03 r_2- The same differences in mass profile shape 
are also easily appreciated in the scaled circular velocity 
profiles, which indicate that the departures from similarity 
are genuine and not just caused by inaccuracies in the scaling 
or by the "bumps and wiggles" caused by unrelaxed tidal 
debris and substructure. 

We have verified this further by performing the same 
analysis after removing bound substructure clumps iden- 
tified by SUBFIND: the same conclusion applies to the 
"cleaned" profiles of the main smooth halo. With hindsight, 
this is perhaps not too surprising. Bound substructures do 
not amount to more than ~ 10% of the halo mass (Springel 
et al., 2008b), and therefore cannot alter the results dis- 
cussed above. 

We have also checked that the differences in a are not 
caused by transient departures from equilibrium or numeri- 
cal resolution: the same qualitative trends, and indeed very 
similar a values, are seen at earlier times and in runs with 
fewer particles. There also seems to be little correlation be- 
tween a and the overall triaxiality of the system; however, 
we shall only deal here with spherically-averaged profiles, 
and defer a detailed study of departures from sphericity to 
a later paper. 

Although the departures from similarity appear signif- 
icant, we must also emphasize that they are rather subtle, 
and are only clearly evident because of the large radial range 
resolved by our simulations, about three decades in radius 
within the virialized region of a halo. Simulations with more 
limited numerical resolution have hinted at this but had dif- 
ficulty making such a compelling case for non-similarity (see, 
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Figure 7. Maximum value of the asymptotic inner slope of the 
density cusp, as a function of radius for our Aq-A convergence se- 
ries. Excellent numerical convergence is achieved at radii compa- 
(7) 

rable to r co ' nv (the inner limit of the thick lines; thin lines extend 
down to r-conv). This shows that there is not enough mass near 
the centre of Aq-A to sustain a cusp steeper than p ex. r~ 0-9±0.l 
Arrows are as in Fig. 1. 

e.g., Navarro et al., 2004; Merritt et al., 2005; Stoehr, 2006; 
Merritt et al., 2006). 



3.6 The Cusp 

It is clear from the residuals in the bottom panels of Fig. 3 
that, near the centre, the M99 profile approximates the sim- 
ulated halos more poorly than either NFW or Einasto. The 
weak performance of the M99 formula may be traced to 
its steep asymptotic inner slope, par" J . Indeed, all six 
Aquarius halos have measured slopes in the inner regions 
that are substantially shallower than —1.5. This is shown 
in Figs. 5 and 6, where the thick portion of each curve cor- 
responds to r > 7"conv an d the innermost point plotted to 
rconv In ah cases, the logarithmic slopes converge well in- 
side rclnv, and only minor deviations may be seen at radii 
beyond r£L- 

Interestingly, the slope of the Aq-A-1 profile at r — rionv 
is exactly —1, and becomes shallower inward, so it is clear 
that at least for this halo we are able to resolve a region 
where the dark matter profile has become shallower than 
— 1, the asymptotic value of the NFW profile. Fig. 6 shows 
the radial dependence of the logarithmic slope for all six 
level-2 halos and confirms the general applicability of the 
Aq-A results: the measured slopes of all halos approach — 1 
(and are certainly shallower than —1.5) at the innermost 
resolved point. 

Figs. 5 and 6 also make clear that there is no sign 
that the profiles are approaching power-law behaviour near 



Figure 8. As Fig. 7, but for our six level-2 Aquarius halos. Re- 
sults are similar in all cases and rule out cusps steeper than r _1 
for ACDM halos. 

the centre: they keep getting shallower to the innermost re- 
solved radius. This behaviour is well captured by the Einasto 
model, where the logarithmic slope is simply a power-law of 
radius, din p/ din r tx r°. Our results thus rule out recent 
claims of cusps as steep as r~ 1,2 in typical ACDM halos 
(Diemand et al., 2004, 2005, 2008). 

This conclusion is unlikely to depend on the details of 
our profile construction and/or fitting procedures. Indeed, as 
we show in the next subsection, there is actually not enough 
mass within the innermost resolved radius to allow for a cusp 
as steep as r -1,2 . Recent work by Stadel et al. (2008), also 
based on very high-resolution simulations, agrees with our 
present conclusions, and argues for asymptotic inner slopes 
shallower than —1, as previously suggested by Navarro et al. 
(2004). 



3.7 The Asymptotic Inner Slope 

The results presented above do not preclude the possibility 
that a shallow power-law cusp may be present in the inner- 
most regions which are still unresolved in our simulations. It 
is therefore interesting to estimate the maximum value that 
the slope of such a cusp may take. This is constrained, at 
any radius, by the total enclosed mass and the local value 
of the spherically averaged density: slopes steeper than 7 max 
require more mass than is available within that radius. This 
constraint assumes only that the logarithmic slope is mono- 
tonic with radius and that the halo is not hollow. It is then 
straightforward to show that the maximum possible inner 
asymptotic slope is 7 max = 3(1 — p(»")/p(r)), where p(r) is 
the mean density enclosed within r. Evaluated at the inner- 
most radius where both local density and enclosed mass (or, 
equivalently, circular velocity) have converged, this quan- 
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tity provides an important constraint on the density profile 
at radii that remain unresolved even in our best simulations. 

We show this parameter as a function of radius for our 
Aq-A convergence series in Fig. 7. This figure shows that 
7max converges to better than 0.1 for r > rioLv (the inner- 
most point of the thick portion of the profiles). Our data for 
Aq-A thus indicates that there is not enough mass in the 
unresolved region to support a cusp steeper than r - - 9 ± 01 . 
Fig. 8 shows that the results for Aq-A are not exceptional: 
all our level-2 Aquarius halos suggest maximum possible 
asymptotic slopes of about —1. 



4 DYNAMICAL PROFILES 

4.1 Velocity Dispersion Structure 

Fig. 9 shows velocity dispersion and anisotropy profiles for 
our Aq-A series and demonstrates that the excellent numer- 
ical convergence of our simulations extends to their veloc- 
ity dispersion structure. The velocity dispersion (squared) 
is computed simply as twice the specific kinetic energy in 
each spherical shell and the anisotropy as /3 = 1 — a 2 /(2a 2 ), 
where of and a 2 are the (squared) velocity dispersion in 
tangential and radial motions, respectively. Besides numer- 
ical convergence, the panels in this figure illustrate two im- 
portant points. The first concerns the shape of the velocity 
dispersion profiles (left panel in Fig. 9), which is remarkably 
similar to that of the r 2 p profiles shown in Fig. 1. This co- 
incidence suggests an intimate connection between density 
and velocity dispersion, which we explore in more detail in 
Sec. 4.4. The second point concerns the anisotropy profile, 
which is clearly non-monotonic. It is nearly isotropic at the 
centre, becomes radially anisotropic at intermediate radii, 
but the dominance of radial motions decreases again near 
the virial radius. As shown in Fig. 10, these properties ap- 
pear to be rather general, since all six Aquarius halos have 
non-monotonic anisotropy profiles and similar velocity dis- 
persion profile shapes. 

4.2 Self-similarity? 

Fig. 10 also demonstrates a clear lack of self-similarity in the 
structure of the simulated halos. We have chosen to empha- 
size this by rescaling all profiles so as to match the peak of 
the a(r) curve, which occurs at r(er max ). This scaling demon- 
strates that, as with the density profiles, the shape of the 
a(r) profiles differs subtly but significantly amongst halos. 
We have checked that these differences in shape are not due 
to bound subhalos; removing all the subhalos identified by 
our SUBFIND algorithm and recalculating the dispersion 
and anisotropy profiles results in only rather minor changes 
The most striking case is again that of halo Aq-E-2 (blue 
curve), whose a(r) profile is much broader than the others. 
Recall that this halo also stands out in Fig. 3 as having an 
unusually broad r 2 p profile. Halo Aq-E-2 also has an unusual 
velocity anisotropy profile, with less predominance of radial 
motions than the rest of the series. The departures from sim- 
ilarity in mass and velocity structure therefore seem closely 
linked, suggesting that these halos may share a common 
property that combines density and velocity dispersion. We 
explore this in Sec. 4.4 below. 



4.3 Anisotropy-slope relation 

We may use the results of the previous subsection to assess 
recent claims by Hansen & Moore (2006) of a general con- 
nection between the local values of logarithmic slope, 7, and 
the velocity anisotropy, f3. We show this in Fig. 11, where 
we plot P vs 7 for all level-2 Aquarius halos. Open circles 
correspond to the inner regions of the halo (r^ v < r < r-2) 
whereas filled circles correspond to the outer regions (r_2 < 
r < r2oo). As in other figures, different colours correspond 
to the different Aquarius halos. The relation proposed by 
Hansen & Moore is shown by a dashed line and accounts 
reasonably well (albeit not perfectly) for our data in the in- 
ner regions where both the anisotropy and the logarithmic 
slope are monotonic functions of r. 

However, there are large departures from this relation in 
the outer regions, where the density profile steepens further 
but the velocity ellipsoid tends to become less anisotropic. 
The failure of the Hansen & Moore relation in the outer 
regions is not unexpected since 7, unlike /3, is monotonic 
with radius. We conclude that, if a simple relation links 
anisotropy and slope, it can only hold in the inner regions 
of halos. 



4.4 The Phase-Space Density Profile 

The similarity in shape between the a 2 and r 2 p profiles high- 
lighted above suggests that there may be a simple scaling 
between densities and velocity dispersions in halos. This is 
best appreciated by considering the quantity p/a 3 , which, 
for dimensional reasons, we shall call the pseudo-phase-spacc 
density, although it is important to realise that it is not 
the true coarse-grained phase-space density at the resolu- 
tion of our simulations, or even the average of this quantity 
in spherical shells. For consistency with the rest of our anal- 
ysis, we calculate p/a 3 directly from the estimates of p and 
a computed in concentric spherical shells. 

Fig. 12 shows the p/a 3 profile for our Aq-A conver- 
gence series. As noted by Taylor & Navarro (2001), the pro- 
file of this quantity is remarkably well approximated by a 
power-law. More remarkable still is the fact that the power 
law is indistinguishable from that predicted by the similar- 
ity solution of Bertschinger (1985) for infall onto a point 
mass in an otherwise unperturbed Einstein-de Sitter uni- 
verse, p/a 3 oc r~ 1,875 (dot-dashed line in Fig. 12). This so- 
lution is spherically symmetric, involves purely radial mo- 
tions, and is violently dynamically unstable, so its relevance 
to ACDM halos is far from clear. The residuals in the bot- 
tom panel of Fig. 12 are deviations from a Bertschinger law 
matched within the characteristic radius r_2, where sub- 
structure bumps and wiggles are minimal. 

Note that, although there is only one free parameter in 
this fit (the vertical scaling), the residuals do not exceed 
~ 20% anywhere within the virial radius, even though sub- 
structures add significant noise to the dynamical measure- 
ments in the outskirts of the halo. Interestingly, the residuals 
increase when ay, the velocity dispersion in radial motions, 
is used in place of the full 3D rms velocity, a, to estimate the 
"phase-space density". Thus, the r -1,875 behaviour seems to 
concern the full kinetic energy content of each shell rather 
than just radial or tangential motions. 

Fig. 13 shows that similar conclusions apply to the rest 
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Figure 9. Left panel: Velocity dispersion profiles for our Aq-A convergence series. Arrows, line-types and colours are as in Fig. 1. Note 
the excellent numerical convergence. The shape of the velocity dispersion profile is remarkably similar to that of the r 2 p profile shown in 
Fig. 1, highlighting the intimate connection between the density and velocity dispersion profiles which is responsible for the power-law 
behaviour of the pseudo-phase-space density profile discussed in Sec. 4.4. Right panel: Anisotropy profiles for the Aq-A convergence 
series. Note the non-monotonic variation with radius: the halo is nearly isotropic near the centre, is dominated by radial motions at 
intermediate radii, but becomes markedly less anisotropic near the virial radius. 




Figure 10. As Fig. 9, but for all six level-2 resolution Aquarius halos, scaled to match at the peak of the profile, identified by <T ma x 
and r(cr max ). This scaling highlights small but significant departures from similarity in the velocity dispersion structure of ACDM halos. 
Note the correspondence in shape between the velocity dispersion and r 2 p profiles shown in Fig. 1, which reflects the "universal" pseudo- 
phase-space density profile of the halos (Fig. 13). Note also that the non-monotonic behaviour of the anisotropy highlighted in Fig. 9 is 
common to all six halos. 
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of the Aquarius halos. Residuals from the Bertschinger law 
are small for all halos, and are typically larger when the 
radial velocity dispersion is used. Note that there is some 
"curvature" in the residual profiles, suggesting that a power- 
law is a good, but perhaps not perfect, description of the 
radial dependence of p/a 3 . We are currently investigating 
the origin of this curvature and plan to report on it in a 
future paper (Ludlow et al., in preparation). 

A power-law radial dependence is approximately pre- 
served when o r is used, but the best fitting value of the 
exponent differs systematically from —1.875. This may be 
seen in the bottom panels of Fig. 13, which show the resid- 
uals from the best fitting p/a 3 oc r x law. The values of the 
best-fit exponent for both p/a 3 and p/a 3 (x and Xr, respec- 
tively) are listed in Table 2. 

Perhaps the most important result from Fig. 13 is that 
there seems to be very little scatter between halos when 
considering their p/a 3 profiles. Take, for example, the case of 
halo Aq-E-2, which was a clear outlier in the density, velocity 
dispersion, and anisotropy profiles. When considering pjo 3 
this halo is unremarkable, and follows the Bertschinger law 
as closely as the others. 

This shows that there is a sense in which ACDM halos 
are nearly universal, but that universality does not extend 
to their density or velocity dispersion profiles separately, but 
rather only to their pseudo-phase-space density profile. This 
may appear a bold statement, and it certainly needs to be 
corroborated by future work, but it offers an intriguing per- 
spective into the origin of the near-universal density profiles 
of halos, the meaning of the Einasto shape parameter, a, 
and the provenance of their velocity dispersion structure. 
These issues deserve further investigation. 

We end by noting that, although it is still not clear 
what leads to the power-law stratification of p/a , these re- 
sults may be used to place constraints on the structure of 
the central cusp, under the plausible (but admittedly un- 
proven) assumption that the power-law behaviour of the 
phase-space density continues all the way to the centre. For 
example, Taylor & Navarro (2001) used this assumption to 
show that, for isotropic systems, a power-law pseudo-phase- 
space density implies an inner density cusp with p oc r -0 ' 75 . 
This is certainly consistent with the results shown in Fig. 7, 
which only exclude cusps steeper than r -° 9 ± 01 However, 
as we show in Fig 5, the detailed profile which they derive for 
an isotropic halo with Bertschinger's power-law p/a 3 profile 
is a significantly worse fit to our numerical data than the 
Einasto profile. 

The power-law behaviour of the pseudo phase-space 
density has been confirmed by a number of authors, and 
seems to be present even at early redshift (Vass et al., 2008). 
Interestingly, the average power-law exponent to the p/a 3 
profile is (x r ) ~ 1.97, close to the "critical" 1.94 required by 
Dehnen & McLaughlin (2005) to have a dynamical model 
that is well behaved at all radii. Simulations of even larger 
dynamic range seem required in order to explore the true 
asymptotic inner behaviour of the dynamical profile of a 
halo, if indeed there is any such asymptote. 
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Figure 11. Local values of the logarithmic slope of the density 
profile plotted versus velocity anisotropy. The relation proposed 
by Hansen & Moore (2006) is shown as a dashed line. Because 
the density profile steepens gradually from the centre outwards 
whereas the velocity anisotropy is non-monotonic, no simple re- 
lation between these two quantities is valid throughout the halos. 
The Hansen & Moore formula approximates our results quite well 
in the inner regions, but large deviations may be seen outside r_2, 
particularly at the largest radii where our halos are approximately 
isotropic but their density profiles are steepest. Open circles cor- 
respond to r^onv < r < r_2, filled circles to r_2 < r < V2oo- 
Colors are as in Fig. 3. 



5 SUMMARY 

We have analysed density, velocity dispersion, anisotropy 
and pseudo-phase-space density profiles at redshift zero for 
simulated halos from the Aquarius Project. This is a set of 
six galaxy-sized halos whose formation and evolution have 
been simulated at a variety of resolutions in their proper 
ACDM context. The set includes the largest simulation of 
this kind reported so far; a ~ 4.4 billion particle simulation 
in which the final halo has 1.1 billion particles within its 
virial radius, r^oo- The set also includes simulations of all 
six halos with 100 - 200 million particles within the virial 
radius, as well as a comprehensive numerical convergence 
study for the largest system. Our main conclusions are as 
follows. 

• Density profiles deviate slightly but significantly from 
the NFW model, and are approximated well by a fitting 
formula where the logarithmic slope is a power-law of ra- 
dius: the Einasto profile (eq. 4). The steeply-cusped profile 
of Moore et al. (1999) is a poor fit to the structure of our 
six halos. 

• We find convincing evidence that the shape parame- 
ter of the Einasto formula varies from halo to halo at given 
mass (see Table 2) . This complements the earlier conclusion 
of Merritt et al. (2006), Gao et al. (2008) and Hayashi & 
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Figure 12. Pseudo-phase-space density profiles for our Aq-A convergence series, estimated as p/o- 3 , computed in concentric spherical 
shells. Arrows, line-types, and colours are as in Fig. 1. Note the remarkable power-law behaviour of this quantity, a result already noted by 
Taylor & Navarro (2001). The dot-dashed line is not a fit to the data, but rather the prediction of the similarity solution of Bertschingcr 
(1985) for infall onto a point mass in an otherwise unperturbed Einstein-de Sitter universe, p/cr 3 ex. r -1 - 875 . This has been scaled to 
match Aq-A at r < r_2- Residuals from the Bertschinger solution are shown in the bottom panels. Note that this power-law behaviour 
is most evident when the full 3D velocity dispersion is used (left panels). When only the radial velocity dispersion is used (right panels) 
deviations from the Bertschinger solution are considerably larger. 



White (2008) that its mean value varies systematically with 
halo mass. Together these results imply that the density 
profiles of ACDM halos are not strictly self-similar: different 
halos cannot be rescaled to look alike. This lack of similar- 
ity extends to the kinematic structure, as measured by the 
velocity dispersion and anisotropy profiles. 

• Intriguingly, departures from similarity are minimized 
when analyzing a pseudo-phase-space density profile defined 
as p/cr 3 . This suggests a limited sense in which ACDM halos 
are indeed nearly "universal" . The pseudo- phase-space den- 
sity profiles are very well approximated by p/cr 3 oc r -1 ' 875 , 
the power law predicted by Bertschinger's similarity solution 
for infall onto a point mass in an otherwise unperturbed 
Einstein-de Sitter universe. This simple law has only one 
scaling parameter and no shape parameters, yet it approxi- 
mates, for over six decades, the p/cr 3 profiles to better than 
20-30%, all the way from the innermost resolved point to 
the virial radius. The power-law description is, however, not 
perfect, and further work designed to understand better its 
origin and limitations seems warranted. 

• Density profiles become monotonically shallower in- 
wards, down to the innermost resolved point, with no indi- 
cation that they approach power-law behaviour. The inner- 
most slope we measure is slightly shallower than —1, a result 
supported by estimates of the maximum possible asymptotic 
inner slope. 

• These results convincingly rule out recent claims that 
typical ACDM halos may have asymptotic central cusps as 
steep as r" 12 (Diemand et al., 2004, 2005, 2008). Shallower 



cusps, such as the asymptotic r behaviour predicted 
by the model of Taylor & Navarro (2001), cannot yet be 
excluded. These results should discourage further work as- 
suming CDM cusps steeper than r _1 except possibly around 
central black holes. 

• Velocity anisotropy does not depend monotonically on 
radius beyond r_2. Halos are roughly isotropic near the cen- 
tre, are dominated by radial motions at intermediate radii, 
but become more isotropic again as the virial radius is ap- 
proached. This behaviour does not appear to be driven by 
the presence of substructure. Given that the slope of the 
density profile does increase monotonically with radius, this 
implies that no simple relation between anisotropy and slope 
can hold throughout a halo. The relation recently proposed 
by Hansen & Moore (2006) works reasonably well in the 
inner regions (r < r-2), but fails at larger radii. 

The main aim of the Aquarius Project is to provide reli- 
able theoretical predictions for the structure and formation 
history of dark matter halos like that surrounding the Milky 
Way down to radii of order 100 pc. This permits direct com- 
parisons with a number of observations with minimal extrap- 
olation, and it helps to design new observational strategies 
aimed at testing the cold dark matter paradigm on these 
very non-linear scales. 

We recognize, however, that many of these tests and 
predictions will apply to regions where baryons play an im- 
portant dynamical role. Our numerical work provides robust 
results for the limiting but unrealistic case of pure dark mat- 
ter halos, and these will undoubtedly be modified in non- 
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Figure 13. Pscudo-phasc-spacc density profiles of all six level-2 Aquarius halos. Radii have been scaled to r_2, and the pseudo-phase- 
space densities to maximise agreement within r_2. Note that for all six halos these profiles are very well approximated by power laws 
with an exponent very close to that of the Bertschinger solution. All halos, including those that were outliers in the density, velocity 
dispersion, and anisotropy profiles, are almost indistinguishable in this plot. Deviations from the Bertschinger law are typically more 
pronounced when radial velocity dispersion is used instead of the full 3D velocity dispersion. Residuals from the best-fit power-laws, 
p/u 3 oc r x , are shown in the bottom panels. The values of \ are listed for each halo in Table 2. 



trivial ways by the presence of baryons. Providing a full 
account of the coupled structure of the cold dark matter 
and baryonic components in galaxies like our own is clearly 
the next major computational challenge, and it is likely to 
exercise us for some time to come. 



ACKNOWLEDGMENTS 

The simulations for the Aquarius Project were carried out at 
the Leibniz Computing Center, Garching, Germany, at the 
Computing Centre of the Max-Planck-Society in Garching, 
at the Institute for Computational Cosmology in Durham, 
and on the 'STELLA' supercomputer of the LOFAR ex- 
periment at the University of Groningen. This work was 
supported in part by an STFC rolling grant to the ICC. 
CSF acknowledges a Royal Society Wolfson Research Merit 
award. AH acknowledges financial support from NOVA and 
NWO. 

REFERENCES 

Allgood B., Flores R. A., Primack J. R., et al., 2006, MNRAS, 

367, 1781 
Bertschinger E., 1985, ApJS, 58, 39 

Dehnen W., McLaughlin D. E., 2005, MNRAS, 363, 1057 
Diemand J., Kuhlen M., Madau P., 2007, ApJ, 657, 262 



Diemand J., Kuhlen M., Madau P., et al., 2008, Nature, 454, 735 
Diemand J., Moore B., Stadel J., 2004, MNRAS, 353, 624 
Diemand J., Zemp M., Moore B., Stadel J., Carollo C. M., 2005, 

MNRAS, 364, 665 
Einasto J., 1965, Trudy Inst. Astroz. Alma-Ata, 51, 87 
Eke V. R., Cole S., Frenk C. S., 1996, MNRAS, 282, 263 
Frenk C. S., White S. D. M., Davis M., Efstathiou G., 1988, ApJ, 

327, 507 

Fukushige T., Makino J., 2001, ApJ, 557, 533 

Gao L., Navarro J. F., Cole S., et al., 2008, MNRAS, 387, 536 

Gao L., White S. D. M., Jenkins A., Stoehr F., Springel V., 2004, 

MNRAS, 355, 819 
Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel 

J., 2000, ApJ, 544, 616 
Graham A. W., Merritt D., Moore B., Diemand J., Terzic B., 

2006, AJ, 132, 2701 
Hansen S. H., Moore B., 2006, New Astronomy, 11, 333 
Hayashi E., Navarro J. F., Springel V., 2007, MNRAS, 377, 50 
Hayashi E., White S. D. M., 2008, MNRAS, 388, 2 
Henry J. P., Evrard A. E., Hoekstra H., Babul A., Mahdavi A., 

2008, ArXiv e-prints 
Jing Y. P., Suto Y., 2002, ApJ, 574, 538 

Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ, 
522, 82 

Komatsu E., Dunkley J., Nolta M. R., et al., 2008, ArXiv e-prints, 
0803.0547 

Kuhlen M., Diemand J., Madau P., 2008, ArXiv e-prints, 
0805.4416 

Merritt D., Graham A. W., Moore B., Diemand J., Terzic B., 



The Diversity and Similarity of Simulated Cold Dark Matter Halos 15 



Halo 


m-p 




^200 


M 20 o 


JV200 


^max 


^ max 


CT host 


*7 max 




[M Q /h] 


[pc/h] 


[kpc/h] 


[M @ /h] 


[106] 


[km/s] 


[kpc/h] 


[km/s] 


[km/s] 


Aq-A-1 


1.250x10 s 


14 


179.41 


1.343X10 12 


1074.06 


208.75 


20.69 


117.47 


261.70 


Aq-A-2 


l.OOOxlO 4 


48 


179.49 


1.345X10 12 


134.47 


208.49 


20.54 


117.13 


261.88 


Aq-A-3 


3.585X10 4 


87 


179.31 


1.341X10 12 


37.39 


209.22 


20.35 


117.31 


262.80 


Aq-A-4 


2.868 x10 s 


250 


179.36 


1.342X10 12 


4.68 


209.24 


20.58 


117.23 


262.29 


Aq-A-5 


2.294X10 6 


500 


180.05 


1.357X10 12 


0.59 


209.17 


20.84 


116.61 


260.59 



Aq-B-2 


4.706 xlO 3 


48 


137.02 


5.982x10" 


127.09 


157.68 


29.31 


89.59 


190.74 


Aq-C-2 


1.021 xlO 4 


48 


177.26 


1.295X10 12 


126.77 


222.40 


23.70 


124.08 


270.50 


Aq-D-2 


1.020X10 4 


48 


177.28 


1.295X10 12 


126.98 


203.20 


39.48 


113.15 


254.28 


Aq-E-2 


7.002 xlO 3 


48 


154.96 


8.652X10 11 


123.56 


179.00 


40.52 


101.73 


215.14 


Aq-F-2 


4.946 xlO 3 


48 


152.72 


8.282X10 11 


167.45 


169.08 


31.15 


96.78 


204.53 



Table 1. Basic parameters of the Aquarius simulations. We have simulated 6 different halos, each at several different numerical resolutions. 
The leftmost column gives the simulation name, encoding the halo (A to F), and the resolution level (1 to 5; 1 is our highest resolution, 
5 the lowest). m p is the particle mass in the high-resolution region, eq is the Plummer-equivalent gravitational softening length, r2oo is 
the virial radius, defined as the radius enclosing a mean overdensity 200 times the critical value for closure, M200 is the mass within the 
virial radius, ./V200 is the total number of particles within r2oo- Other characteristic properties of the halos listed are the position (r max ) 
of the peak (V m ax) of the circular velocity profile, as well as the ID velocity dispersion of the main halo (cr host ), and the peak (<r m ax) of 
the velocity dispersion profile. 
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Table 2. Fit parameters of Aquarius halos. The first column labels each halo, as in Table 1, the second and third list the convergence 
radii obtained for k = 1 and k = 7. These radii, r'onv an d r£j, v , respectively, correspond to where departures from convergence in the 
circular velocity are expected to be of order 10% and 2.5%. The characteristic scale radius r_2 corresponds to where the logarithmic 
slope equals the isothermal value; p_2 = p{r-2), and a is the best-fit Einasto parameter. The uncertainty in a indicates the range where 
AQ/Q deviates by less than 50% from the absolute minimum shown in Fig. 4. Strictly, these are non-symmetric, so we conservatively 
quote the largest deviation, positive or negative. \ refers to the exponent of the best fitting power-law describing the p/cr 3 profile. \r 
is analogous to x> but for p/<7 3 , where a r is the rms velocity in radial motions. \ and \ r are computed by minimizing residuals in the 
region r^ nv < r < r_2- Finally, 7 max lists the value of the maximum asymptotic slope of the density profile cusp, measured at r = r^} nv . 



2006, AJ, 132, 2685 
Merritt D., Navarro J. F., Ludlow A., Jenkins A., 2005, ApJL, 
624, L85 

Moore B., Ghigna S., Governato F., et al., 1999a, ApJL, 524, L19 
Moore B., Quinn T., Governato F., Stadel J., Lake G., 1999b, 

MNRAS, 310, 1147 
Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563 
Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493 
Navarro J. F., Hayashi E., Power C, ct al., 2004, MNRAS, 349, 

1039 

Power C, Navarro J. F., Jenkins A., et al., 2003, MNRAS, 338, 
14 

Prugnicl P., Simien F., 1997, A&A, 321, 111 

Sersic J. L., 1968, Atlas de galaxias australcs, Cordoba, Ar- 
gentina: Observatorio Astronomico, 1968 
Spergel D. N., Verde L., Peiris H. V, et al., 2003, ApJS, 148, 175 



Springel V., 2005, MNRAS, 364, 1105 

Springel V., White S. D. M., Jenkins A., et al., 2005, Nature, 435, 
629 

Springel V., Yoshida N., White S. D. M., 2001, New Astronomy, 
6, 79 

Springel et al., 2008a, ArXiv e-prints, 0809.0894 

Springel et al., 2008b, ArXiv c-prints, 0809.0898 

Stadel J., Potter D., Moore B., et al., 2008, ArXiv e-prints, 808 

Stoehr F., 2006, MNRAS, 365, 147 

Stoehr F., White S. D. M., Springel V., Tormen G., Yoshida N., 

2003, MNRAS, 345, 1313 
Taylor J. E., Navarro J. F., 2001, ApJ, 563, 483 
Vass I., Valluri M., Kravtsov A., Kazantzidis S., 2008, ArXiv c- 

prints 

Vogelsberger M., Helmi A., Springel V., et al., 2009, MNRAS, 
395, 797 



16 Navarro et al. 



White S. D. M., 1996, in Cosmology and Largc-Scale Structure, 
edited by R. Schaefer, J. Silk, M. Spiro, J. Zinn-Justin, Dor- 
drecht: Elsevier, astro-ph/9410043 



